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Abstract 

A method for the fast and accurate solution of the radiative trans- 
fer equation in plane-parallel media with coherent isotropic scattering is 
presented. This largely analytical approach uses the formalism of mero- 
morphic functions in order to obtain the total solution, i.e. the angle and 
depth variation of the radiation field. A discretization of space and angle 
coordinates is not required. For the further application in the accretion 
disks a particular case of a finite slab whose the midplane is the symmetry 
plane is considered. 



1 Introduction 

The transport of energy by means of radiation plays a very important role 
in, e.g., plasma physics, environmental physics, and especially in astrophysics. 
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The quantity most commonly used to describe a radiation field is the specific 
intensity, which obeys the radiative transfer equation (RTE) (cf. [Q). In most 
cases this is an intcgro-differential equation. In many applications (e.g. stars) 
it is sufficient to consider a ID radiative transfer equation. The plane-parallel 
RTE can also be useful in modeling of accretion disks. In general, an accretion 
disk is a complicated 3D structure and should be treated as such, however the 
multi-dimensional RTE remains computationally costly. An alternative method 
used by various authors [| |, | is to represent an accretion disk as a system 
of independent rings, each of them radiating as a plane-parallel slab. Although 
such an assumption is unrealistic, it nevertheless represents a landmark test in 
the theory of accretion disks. 

The ID radiative transfer equation is rather well developed and there exists a 
vast literature on this topic ( see e.g. (|, @, |[ |§1 ) ■ Among numerical methods 
the discrete ordinate method is the most usable. The approximation of the 
definite integral in the RTE by a quadrature sum leads to the replacement 
of the continuous radiation field by a finite set of pencils of radiation. The 
subsequent discretization of the transport operator results in a system of linear 
equations which is usually solved itcratively by Jacobi iteration ("accelerated 
A-iteration"). This accurate and fast method fails, however, in media with large 
optical depths and/or a large scattering fraction. Problems may also appear in 
media with steep gradients. 

Methods which proceed analytically as far as possible are most desirable. 
Usually they give more insight into the general behavior of the solution and are 
less CPU-time and memory consuming. A method developed by Chandrasekhar 
[|| can be included in this class of methods. It was found that the exact solu- 
tion of the RTE in semi-infinite atmospheres leads to closed expressions for the 
angular distribution of the emergent radiation, involving a so-called iJ-function 
that is the solution of an integral equation of standard form. In media of fi- 
nite optical depth the emergent radiation can be expressed in terms of certain 
rational functions X and Y which satisfy integral equations too. The method 
is not commonly used because the Chandrasekhar functions H , X and Y are 
difficult to calculate in spite of the different methods proposed [|o[ 0, |l2|, |l^ . 
In addition, it does not provide the distribution of the specific intensity with 
depth that is required in various models of astrophysical objects. 

In the present paper we continue discussing a method for the analytical 
solution of the plane-parallel RTE in finite media. In the previous papers ( |p^ , 
|l5| ) the attention was mostly paid to mathematical aspects of the problem and 
algorithmic aspects were hardly considered. Only the angular distribution of the 
emergent intensity was found and only in the lowest " separable" approximation 
(see below). In this paper we solve an extended problem: we also look for the 
solution of the RTE at every point in a slab. We have slightly altered initial 
conditions and now consider a medium with no incident radiation but with a 
non-zero depth-dependent source function. The midplane of the medium is the 
symmetry plane which implies a further application of the solution to accretion 
disks. 

By introducing intensities in outward and inward directions we are able to 
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write the RTE in a matrix form. The involved matrix M is an operator in 
an infinite dimensional space and possesses rather unpleasant properties, in 
particular, its spectrum extends from —oo to oo. In spite of this, the formal 
solution can be expressed in terms of hyperbolic functions of the matrix M. 
Note that these functions are bounded and therefore no difficulties arise related 
to the infinite spectrum of M . The successive application of the formalism of 
meromorphic functions and Krein's formula for the finding of inverse operators 
leads to the exact solution represented in a form of infinite series. Because 
of their very slow convergence these results are not well suited for numerical 
work. Therefore an appropriate approximation of the infinite sums by finite 
ones is proposed where the number of terms defines the order of the separable 
approximation. Experience has shown that the final results obtained in such a 
matter in the lowest approximation are quite reasonable for practical purposes. 
In the 5-6th approximation orders the results can be regarded as highly precise. 



2 Formulation of the problem 

The specific intensity X(t, /i) in plane-parallel media with coherent isotropic 
scattering is described by the following radiative transfer equation (cf. |^]) 

^dX^ = -I(r, m) + /3 I(r, /i') d^^' + eB{T) (1) 

where B(t) is the Planck function, ji = cos0, 6 is the angle between the normal 
and the beam of radiation, r is optical depth, /3 = (1 — e)/2. The de-excitation 
coefficient e can vary in the interval [0, 1]. We assume it is constant. 

Our aim is to find the solution of Eq. (]^), i.e. both the depth and angle 
variations of the radiation field in a slab of finite optical depth. We assume 
that the symmetry plane is at r = 0. The optical depth is measured away from 
the symmetry plane and equals A and —A at the upper and lower boundaries, 
respectively. 

For the boundary conditions we require that no radiation is incident on the 
slab from outside, i.e. 

J+(A,Ai)=J-(-A,/i) = (2) 

where I"*" (r, /z) and I~ (r, /i) are intensities in the positive and negative direc- 
tions with respect to /i, i.e. 

^^{'T,!-!') =I{T,fi), X"(t,^) =X(t,-^) for /i > 0. 



3 Matrix form of the solution 

It is convenient to solve this problem in the matrix formalism. The representa- 
tion of the specific intensity as a two-component vector 
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enables us to write Eq. (|I|) in compact matrix notation 
dl(r, fi) 



dT 



-M^^.I(t,a/)+B(t) (3) 



where the integration over the index fi' is implied. The matrix M and the vector 
B are given by 

M^^, = B^,, - (|v) (3 (u|)^^,, B(r) = eB(r)|v) 

with 

Da'm' = ^-^s^Cm- |v) = --^|e_), |u) = -^|e+), 



T3 



The matrices M and T3 do not commute but satisfy the condition 

TsM^ = Mt3 

Subsequently we shall drop the index /i in the notation and use the following 
abbreviation 

FiM)W) = f ^F{M),,,\e.), 
Jo yfJ- 

(u|F(M)|v) = //'^(e+|F(MW|e->- 
J Jo Vf^f^ 

The boundary conditions can be represented in the following form 

- ( '0" ) ' i*-^) - ( L 



or 



1(A) = /out ^^|e+), I(_A)=/out^^|e+) (4) 

where /out = /out(M) is the outgoing intensity. 

The formal solution of Eq. can be written in two equivalent forms 

r 

I(r) = e-^("+^)l(-A) + J e-^^^-^"> B{t') dr' (5) 

-A 
A 

= e-^("-^)l(A) - J e-^(---') B(t') dr' 



which resuhs in 

A 

1(A) = e-2MAj(_^)^ j e-M('^-')B(r')dr'. (6) 

-A 

Muhiplying both sides of (||) by e'^'^, substituting (|) and using B{t) = B{-t) 
we obtain 



The internal distribution of the radiation field is given by the solution of Eq. 
(P. It is sufficient to solve this equation for X+, since due to the symmetry X~ 
follows immediately 

/+(-r)=J-(r). (8) 
Thus for the given mean intensity and the boundary conditions we have 

r 

X+(r, /.) = / [(1 - e) J(r') + eB(r')] — (9) 

-A 

The mean intensity is given by 
1 



j{t) = -y^ {T+{T^^,')+l-{T,^l'))d^i' 

= \ f ^{I+{r,^,') + ^{T,^,'))d^,' ^-{n\{i+{r) + r{T))\e+). 

^ Jo vA* ^ 
Using the definition of I(t) and Eq. we obtain 

I(r)+I(-r) = (/+(r)+/-(r))|e+). 

So that the mean intensity becomes 

J(r) = i(u|I(r)+I(-r)). (10) 
Using Eq. (||) and the relation 

1(A) +I(-A)= /out |e+) 

we obtain 

-Mr 



f I -M(-A+r-r') -M(A+r-r') | 

/ f - 2co,h(MA) - - -' 2co,h(MA) ^f^'' * 



-A 



5 



where the unit step function Q{x) is given by 

e(x) 

The substitution of Eqs. (|ll|) and ^ into leads to 



for a; < 0, 

1 for X > 0. 



where 



J(r) = 

Gi{T,r') = 
G2{t-t')-- 



(Gi(r,T')+G2(T-r')) i3(T')dr' 



(12) 



cosh(MT) 



cosh(MT') 



u 



cosh(MA) Ta + tanh(MA) cosh(MA) 
sinh(M(A- \t-t'\)) 



cosh(MA) 

4 Evaluation of the matrix elements 

In order to evaluate (^ and ( [l^ ) we apply the formalism of mcromorphic func- 
tions and Krein's formula. 

By definition, a meromorphic function is a function that is analytic, except 
for a set of poles For a meromorphic function F{z) decreasing for \z\ — > oo 
the following representation is valid 



Pi-) = E 



Frr 



(13) 



where 



lim [z - £_„i) F{z). 



In our case we have the meromorphic functions 
cosh(]V[a) sinh(Ma) 



a| < A 



cosh(MA) ' cosh(MA) ' 
which in accordance with (^) can be represented by 

p °° P 



F(M) 



E 



iym - |v) /3 (u| 



(14) 



where the constants F^^ depend on the particular form of the function F(M), 
and the simple poles are 



6, 



in f 1 



A \2 
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In order to calculate an inverse operator such as that in ( p^ we use Krein's 
formula which states that for an operator S acting in an appropriate space C 
and being of the following form 



N 



S = H-J2 mc^,{W,\ =H- \V)c{W\ 



the inverse operator S ^ can be represented by 



(15) 



where 



V, 



H is an operator, \Vi) and \Wi) are vectors in the space C, Cij is a number 
matrix. 

The application of Krein's formula gives 



Cm = 1 - ( u 

»i 



l+Ai^^m 



FiM) = F{B) + J2^ 

where 

C(iy„i) 

= 1-2/3 
Using (|6|) and ([l7| ) we obtain 

FiM)\v) 
(u|F(M) 
(u|F(M)|v) 



13F„ 



C{iy„ 



u 



1 



D + iy„ 



D + iy, 
1 



(16) 



(17) 



dfi^ 1-2(3 



arctan(?;„) 



\ ^ Fm 1 

EFm / 



1 



D + iy„ 



In order to evaluate these expressions we apply the following representation 



1 ^ 2to 



1 



where 



p{t) = 



' Pit) 
Ci tt + yl, ' 7o 1 + vU 



W 

(l + /3Hn(i^))' + (^/3i) = 



■dt 



(18) 
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to is a positive root of the equation 



C(to) = l + ^ln 

to 



and 



Since 



dC{t) 



dt 



I -to 
l + to 



to{l-tl) ■ 



E 



F„ 



1 



E 



1 



2/i 

-D2 2fi \n-L> ' fi + n 



we obtain 



u 



sinh(M(A- |T-r'|)) 



cosh(MA) 

cosh(MT) 



cosh(MA) 
cosh(Mr) 



u 



= ^i\r-r'\), 
= 3'(^,/")|v), 



(19) 



where 



and 



$(t,//) 



cosh(MA) 

2 sinh(to(A-T)) 1 /"^ p(t)sinh((A-r)/t) 



pC'i cosh(^oA) f3 Jo t cosh(A/t) 
cosli(T//i) 2to / cosli(T//i) cosli(ioT) 



dt (20) 



cosh(A//i) Ci Vcosh(A//^) cosh(toA) J tg/z^ - 1 

^ 1 / cosh(r/M) _ cosh(r/t) \ m'p(0 
Vcosli(A//i) cosli(A/t)/ _ ^2 



(21) 



The last operator which has to be represented in a form appropriate for 
numerical work is (1 + T3 tanh(MA))^^. The formula ( |l6|) gives the following 
representation (see H)) 



1 + T3 tanh(MA) = + w^^'P^ - u;(2)p+ 

^ (w;(") + )P_ + («;(") - 

where P± = i |e±)(e±| are the projection operators with the properties: 



(22) 



P+ + P- = 1, Pi = P±, p±Pt = 0. 
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The operators w^^^ = w^^^n,fi'), {j = 0, 1,2) are given by 



1 


-I- tanh ^ 


-)1 


- /i') 




4/3 


ni—0 








A 






4(3 


m— 




^/3/2 


A 
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(23) 
(24) 

(25) 



The operators w^^^ and w^"^^ are symmetric, positive definite and have finite 
traces |l^, 15 , i.e. 



Tr w 



(p, fi) d^j. < OO, (j = l,2). 



The inversion of (|2^) is given by 
1 



1 1 
1 + T3 tanh(MA) ^ lyW + iy(o) _ ^^(2) 



so that 



T3 + tanh(MA) 



1 



1 1 



y(0) y;(0)^_^(l) 



(26) 



,(0) 



T3 



5 Separable representation of the solution 



Krein's formula (^_5|) can now be apphed to the expressions in brackets in (^ 
yj{0) corresponds to H and the terms 

and 



^ + y?ntJ''^ l + y2^/i'2' 1 + 1 + y^Ai'^ 

which present in the definitions (|2j) and (^5|) can be regarded as the vectors Vm 
and Wm- Unfortunately, the infinite sums converge very slowly and one needs 
up to several thousands terms in order to achieve the required accuracy. Conse- 
quently the dimension of the matrices involved in (jl^) becomes very large that 
makes the application of Krein's formula inefficient here. In order to accelerate 
computations we use the following approximations for w^^^ and w^^\ 
The operators w*-^-* and w^^) ^^e represented by 



= (27) 



/i^ — fl'- 

■,(2)^, ,,'^1 ^ „3/2 M'" „/3/2 

/9 9 
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where 



A ^„ Cm ' 1 



ni—0 ^ fnr 



(29) 



Let us approximate the function E{t) by a finite sum 



N 



E{t) « E^ii) = ^ 

n=l 

Then we get for <M) and (|8|) 



(30) 



w'^'\^^,^^') « ^K.^i^(M)a„<5„„,V;«(/.') = l^<^V^^'(^'^'l, (31) 

TV 

w'~'\^^,^^') « ^yi2)(/i)a„^„<5„„,F„(^'(M') = |V^^^))j(^^(V^^^)| (32) 



where 



1 + An^i^ ' 



3/2 



- a 5 , 



Z*^' - « /I ^ , 



1 + ' 

The number in (|3^) defines the A^-th separable approximation. 

Taking into aecount the representations ( pl| ) and (32) and using Krein's 
formula (HS) we get now 



1 



1 



1 



1 1 1 



iy(0) _ ^(2) ^{0) ^{Q) 



yd) ^^(1)^^(1) 



1 



1 



w 



(0) 



where 



c(l) ^ 
nn' 



1 7(l)^ c(2) = ^ /(2)^ 



and 



(1) 



^ nrt ' 



1 



2 J (l + A„^2)(l + ^^,^2) 



1 



2y (i + ^„//2)(i + ^^^,^2) 



10 



The substitution of the third term of (pO) into (12) gives 



cosh(MT') 1 cosh(MT) 

T3' 



cosh(MA) w(o) cosh(MA) 



2A dfl 



$(T,/i)$(r',/i)(l + e — ) 



The contribution of the second term of ( pq ) equals zero because P_T3|e_) = 0. 
The substitution of the first term leads to 



u 



cosh(M7 



1 



1 



cosh(MA) \w(°) ~wC^ 



,(0) 



cosh(MT') 



u 



cosh(MT) 



where 



and 



cosh(MA) 



^+'^^cosh(MA) 
cosh(Mr' 



cosh(MA) 



Vr 



(2) 



,(0) 



cosh(Mr) 



cosh(MA) 



cosh(Mr) 



cosh(MA) 



$(t,/i) (1 + e ^ 



1 + Ar^p? 

Taking into account the above expressions the mean intensity becomes 

J(r) = |y'%(|r-r'|)i?(r')dr' + |x„(r)55^\v(r')i?(r')dr' 





2A \ 






he f j 


L 









1 

— < 

and the emergent intensity takes the form 



(33) 



2'out(M) = £(l + e 



- / $(T',^)B(T')dT' 

Jo 



2(1 + A 

where the Einstein's convention is implied. 



K^'{r')B{T')dA. (34) 



6 Discussions 

The practical efficiency of the method depends largely on how well the lower 
orders of the approximation represent the exact solution. In order to obtain 
an accurate solution in the low approximation orders a clever approximation of 
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A=0.1, E=0.5 

0.45 




0.05 I ' ' ' ' ' ' ' ' ' 1 

O 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 

Figure 1: Angular distribution of the outgoing intensity calculated in the dif- 
ferent approximation orders by means of equation with the Planck function 
B{t) = 1 — 0.5 (r/A)^. The curves representing the results of the different 
approximation orders coincide. 



the function E{t) is necessary. The representation of the original function ( p9[ ) 
by the finite sum ( |30| ) with the minimal number of terms and without loss of 
accuracy is the crucial point of our investigation. 

We propose two methods for the approximation (|3^). The first was origi- 
nally developed by Stieltjes and Markov (Appendix A). Using the theory of the 
orthonormal polynomials, the coefficients a„ and An can be found very quickly. 
However, roundoff errors appearing in the calculation of the coefficients pki in 
(A.5) do not allow the method to be implemented, in particular, by means of 
FORTRAN codes. The method is well suited for programs like MATHEMAT- 
ICA. The convergence of the final solution is moderate in optically thick media 
and high in optically thin ones. 

The second method is the so called "Points method". It consist of the 
solution of the system of 2 A'' algebraic equations as described in Appendix B. In 
contrast to the Stieltjes-Markov method it does not suffer from roundoff errors, 
and FORTRAN codes run very well. Some matrices become badly conditioned 
in the high approximation orders {N > 6). However these cases imply a very 
high accuracy of the final solution, which is usually not needed in applications, 
and are not considered. In optically thin media the solution obtained with such 
a„ and An shows approximately the same convergence as in the Stieltjes-Markov 
method. The convergence of the "Points method" is much better in optically 
thick slabs. The calculation of a„ and An takes much longer because one needs 
the values of the original E{t) in the reference points {ti, t2N}- 

Further, all the calculations are carried out with the application of the 
"Points method" where the points are chosen as U = i/2N. 
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A=1000, e=0.5 



0.45 
0.44 
0.43 
0.42 
0.41 
0.4 
0.39 
0.38 
0.37 
0.36 
0.35 




0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 



Figure 2: Angular distribution of the outgoing intensity calculated in the dif- 
ferent approximation orders by means of equation with the Planck function 
B{t) = 1 -0.5(r/A)^ 



Figs. 1^ and ^ show the angular distribution of the emergent intensity in 
slabs of different optical depth. The results were obtained through Eq. (|^) with 
the mean intensity J(t) calculated by means of (js^). The dependence of the 
approximation order is shown. As one can see in optically thin media the precise 
results can be obtained already in the lowest separable approximation whereas 
in optically thick slabs a few additional approximation orders arc necessary in 
order to match the exact solution. 

However, the simplest and the fastest way for the calculation of the emer- 
gent intensity is the direct application of (|3^). As in the previous case a few 
approximation orders are sufRcient to obtain results with reasonable accuracy. 
The exception is for the range of small fi (see Fig. H). The difference between 
this solution and the exact one is caused by errors in the approximation of the 
function E{^^). So, for /i — > we have E{fi^) ^ n, whereas En{ii^) ~ /x^. 
Although the asymptotic behavior of ui*-^-* (/x, and ^^'(/i,^') at small /i are 
the same, w'^^ (/i, /i') and wj^"* (/x, /i') behave in different ways: when ^ 
^(2) ^ while w^^'' ifi, m') ~ V'/'. 

Better results in the region of small fj, can be obtained by means of a further 
improvement of the approximation for the operator w^^'' : 

- VJ^ 

= w^'Hp,^^')+w^'\li,ti') (35) 
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A=10, 6=0.5 



0.55 
0.5 

0.45 
0.4 

0.35 
0.3 

0.25 



O 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 



Figure 3: Angular distribution of the outgoing intensity calculated in the 6-th 
approximation order with the Planck function B{t) = 1 — 0.5 (r/A)^. The 
solid curve represents the precise results obtained by means of (^. Equation 
(^H) was used to get the dotted curve. The results obtained with the improved 
approximation of the operator w*^^^ (|3^) are represented by the dashed-dotted 
curve. 



where 



Applying the "Points method" for the following approximation 

N 



2^ 1^ J_ M _ 



bn 



we get 



AT 



W 



(3) 



St 



hn.B„ 



1 + Bn^l' 



(36) 



(37) 



Thus the final results calculated with the new representation of ui'^^ ( |35| ) show 
better agreement with the exact solution for all ^, as shown in Fig. ^. 

Fig. ^ shows the CPU-time required for the calculation of the internal distri- 
bution of the mean intensity [left panel] Eq . (p3[ )) and the angular distribution 
of the outgoing intensity {right panel] Eq. (|34|)) in the different approximation 
orders. FORTRAN double-precision, optimized codes were used on a HP C240 
computer. Note that time necessary for the calculations of the coefficients a„, 
An and matrix Sl^', is also included. Although these calculations are carried 
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m lout(i^) 




23456 23456 



N-th approximation N-th approximation 

Figure 4: CPU-time required for the calculation of the mean intensity and the 
intensity emerging from the surface of a slab with A = 100 and e = 0.5 in the 
different approximation orders. 



out only once for the given A and e, the most time (about 8 ms in the 2nd 
order) is spent on this, which remains the place for further improvements. 

The comparison of our results with those obtained by the finite element and 
finite difference methods in the case of isothermal media can be found in |16|. 



7 Summary 

In this paper we have given a new method for the solution of the plane-parallel 
radiative transfer equation in media of finite optical depth. The continuous angle 
and depth variation of the specific intensity including the angular distribution 
of the emergent radiation can be obtained fast and accurately for a wide range 
of parameters. The considered configuration of the medium makes the accretion 
disk the most appropriate object for further application of this solution. 

Generalization to the cases of anisotropic and non-coherent scattering is also 
possible. In particular the modeling of a spectral line with complete redistribu- 
tion will be presented in a separate paper. 

We would like to thank K. Manson for valuable comments which have lead 
to an improvement of this paper. This work was supported by the Deutsche 
Forschungsgemeinschaft (Sonderforschungsbereich 359/C2) and the Graduier- 
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tenkoUeg at the Interdisciplinary Center for Scientific Computing of Heidelberg 
University. 



A Approximation of E'(t): Stieltjes-Markov method 

In order to represent E{t) in the form of a finite sum 



ant 



n— 1 

one can apply the Stieltjes-Markov method (see p7[|). 

Let us introduce a measure (f>{dx) on interval [0, 1] in such a way that 

/ ^idx)fix) ^ I -1^ • / f — ^) . (A.2) 

For the approximate evaluation of the integral we use a quadrature formula 
} 

/ ~^c„/(x„) (A.3) 

where Xn are zeros of polynomials of degree A'^ which are orthogonal with respect 
to 4>(dx), i.e. 

1 

(j3{dx) Pm{x)Pi{x) = for m^l, 





and the weights c„ are given by 

^Ar(a;n) 

where Qk{x) are associated polynomials of degree (A^ — 1) 

Qk[s) = / (p{dx) . 

J s-x 



The polynomials Pk (x) can be constructed in the following way. Let us calculate 
2N + 1 moments 



OO / s 1 + fc 



;^ m=0 ^ 



1 / a ^^ tanh (A V^) 
16 



(A.4) 

v=l 



A matrix H with elements 

Hij = ^ 0, ...,7V) 

is symmetric and positive definite. The Cholesky decomposition gives an upper- 
triangular matrix r with the property that H can be written as H = r^r. Then 
the elements 



define the coefficients of the polynomial Pk{x) 

k 



(A.5) 



1=0 

The associated polynomial can be written as 

k-l k 

Qk{x) =^^x'' ^ pkiirii-i-i. 

1=0 i=l+l 

Let us introduce the notation 

s t 



t ■ 



1 + 

1 

1+7 



i-t' 



< s < oo, 



< X < 1, y = 



1-x 



Then we have 



4/3 ^ 1 t _ 4/3 ^ 1 



^ + s 



According to ( |A.2| ) and (A. 3) we get 

1 



E{t) 



(j){dx)- 







C i 



l~x \ X + S 



N 

E 



CnP 



c 



S ~h Xf} 



N 

E 



c„/3 



C ^i^/ i-^" j 2;„ + (1 - a;„)t' 



Thus the coefficients in (A.l) become 

_ c„ 2/3 



A„ = 1. 



(A.6) 
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B Approximation of E{t): "Points method" 

The coefficient a„ and A„ of the A^-th separable approximation can be found 
by solving the system of algebraic equations: 



N 



EN{tl) 



El 



, = ..,27V) 



(B.l) 



where {ti, ...,t2Ar} are points in the interval [0, 1]. The solution of (BJ_) can be 
obtained in the following way. Wc have 



Eiti) 
ti 



N 



N 



k=l 

Let us introduce the notation 

N N 

fc=l s=0 
Mo = 1, U 



Y[{l + Akti) = ^anY[{l + Akti). 

n=l k^n 



■A, 



^ ^ Ai-^ ■ Ai 

l<ii<...<is<N 

N N 

n (1 + Akti) = i?u,u„=o = E ^r'^s 

k^n s— s— 1 

N N 1. 



N 



s-l,,(n) 



AT 



n=l k^n 



Then we can write (B.2) as the following 



AT 



N 



E{ti) 

ti 



s=l s=l 

Let us introduce N x N- matrices and Ai"-component vectors 

f(l) _ f(l) _ -^(^') f(2) _ f(2) _ E{tN+l) 

J ~ Ji ~ J. ^ J h ^ -L 

U tN+l 



(B.2) 



(B.3) 



So that the system of linear equations (B.3) can be rewritten in a form 

Q(i)6-i/(i)u = /W, (B.4) 
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Eliminating b wc obtain the vector u = {ui, ...,U]\[) 

u = [(gW)-ii7(i) - (Q(2))-i^(2) 

The roots of the following polynomial 

N 

i-xf + u,{-xf-^ + ... + itAT = (-1)^ \{{x^ A„) = 

ri=l 

correspond to the constants Ai, An ■ The last step is the solution of the 
matrix equation 

with 

Cis = —^, {l,s^l,...,N) 
t + Asti 

in order to get the coefficients a ~ (ai, cat). 
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